#include "borwein.h"

/* Valores iniciais do método Borwein */
void bo_init_values(mpf_t y, mpf_t a)
{
    /* Calcula a0 = 6-4*sqrt(2) */
    mpf_init_set_ui(a, 2); mpf_sqrt(a, a); mpf_mul_ui(a, a, 4); mpf_ui_sub(a, 6, a);
    
    /* Calcula y0 = sqrt(2)-1 */
    mpf_init_set_ui(y, 2); mpf_sqrt(y, y); mpf_sub_ui(y, y, 1);
}

/* Inicializa estrutura paralela */
void bo_init_param(bo_param *p)
{
    /* Inicia com os valores iniciais */
    bo_init_values(p->y[0], p->a[0]);

    mpf_inits(p->y[1], p->y[2], p->a[1], p->pi[0], p->pi[1], NULL);

    p->stop = B_FALSE;
}

/* Limpa estrutura paralela */
void bo_clear_param(bo_param *p)
{
    mpf_clears(p->y[0], p->y[1], p->y[2], p->a[0], p->a[1], p->pi[0], p->pi[1], NULL);
}

/* Realiza o cálculo de y */
void bo_y_value(mpf_t y1, mpf_t y0)
{
    mpf_t aux1, aux2;
    mpf_inits(aux1, aux2, NULL);

    /* Calcula aux1 = (1-y0^4)^(1/4) */
    mpf_pow_ui(aux1, y0, 4); mpf_ui_sub(aux1, 1, aux1); mpf_sqrt(aux1, aux1); mpf_sqrt(aux1, aux1);
    /* Calcula aux2 = aux1+1 */
    mpf_add_ui(aux2, aux1, 1);
    /* Calcula y1 = (1-aux1)/aux2 */
    mpf_ui_sub(y1, 1, aux1); mpf_div(y1, y1, aux2);

    mpf_clears(aux1, aux2, NULL);
}

/* Realiza o cálculo de a */
void bo_a_value(mpf_t a1, mpf_t a0, mpf_t y1, int k)
{
    int n = 2*k+3;
    mpf_t aux1, aux2;
    mpf_inits(aux1, aux2, NULL);

    /* Calcula aux1 = y1*(1+y1+y1^2) */
    mpf_pow_ui(aux1, y1, 2); mpf_add(aux1, aux1, y1); mpf_add_ui(aux1, aux1, 1); mpf_mul(aux1, aux1, y1);
    /* Calcula aux2 = aux1*2^n */
    mpf_set_ui(aux2, 2); mpf_pow_ui(aux2, aux2, n); mpf_mul(aux2, aux2, aux1);
    /* Calcula a1 = a0*(1+y1)^4 - aux2 */
    mpf_add_ui(a1, y1, 1); mpf_pow_ui(a1, a1, 4); mpf_mul(a1, a1, a0); mpf_sub(a1, a1, aux2);

    mpf_clears(aux1, aux2, NULL);
}

/* Realiza o cálculo de Pi */
void bo_pi_value(mpf_t pi, mpf_t a)
{
    /* Calcula pi = 1/a */
    mpf_ui_div(pi, 1, a);
}

/* Método sequencial Borwein */
void borwein(mpf_t pi, long *it)
{
    mpf_t a[2], y[2];
    mpf_t old_pi;
    int i, i0, i1;

    /* Valores iniciais e inicialização de variáveis */
    bo_init_values(y[0], a[0]);
    mpf_inits(a[1], y[1], old_pi, NULL);

    /* Iterações do algoritmo */
    for(i = 0; ; i++)
    {
        i0 = INDEX2(i, 0);
        i1 = INDEX2(i, 1);

        /* Calcula os parâmetros */
        bo_y_value(y[i1], y[i0]);        
        bo_a_value(a[i1], a[i0], y[i1], i);

        mpf_set(old_pi, pi);
        bo_pi_value(pi, a[i1]);

        /* Verifica convergência */
        if(has_converged(old_pi, pi)) break;
    }

    mpf_clears(a[0], a[1], y[0], y[1], old_pi, NULL);

    *it = i;
}

/* Método paralelo Borwein */
void borwein_conc(mpf_t pi, long *it)
{
    pthread_t thread[3];
    bo_param p;
    int i, j;

    bo_init_param(&p);

    for(i = 0; !p.stop; i++)
    {
        p.i = i;

        /* Cria uma thread para cada parâmetro */
        pthread_create(&thread[0], NULL, bot_calculate_y, &p);
        pthread_create(&thread[1], NULL, bot_calculate_a, &p);
        pthread_create(&thread[2], NULL, bot_calculate_pi, &p);

        /* Espera todas as threads acabar para sincronizá-las */
        for(j = 0; j < 3; j++)
            pthread_join(thread[j], NULL);
    }

    *it = i; /* Número de iterações */
    mpf_set(pi, p.pi[INDEX2(i, 1)]);
    bo_clear_param(&p);
}

/* Borwein Thread para cálculo de y */
void *bot_calculate_y(void *arg)
{
    bo_param *p = (bo_param *) arg;
    int i0 = INDEX3(p->i, 0);
    int i1 = INDEX3(p->i, 1);

    bo_y_value(p->y[i1], p->y[i0]);

    return NULL;
}

/* Borwein Thread para cálculo de a */
void *bot_calculate_a(void *arg)
{
    bo_param *p = (bo_param *) arg;
    int i1, j0, j1;

    /* "Bolhas" do método pipelined */
    if(p->i < 1) return NULL;

    i1 = INDEX3(p->i-1, 1);
    j0 = INDEX2(p->i-1, 0);
    j1 = INDEX2(p->i-1, 1);

    bo_a_value(p->a[j1], p->a[j0], p->y[i1], p->i-1);

    return NULL;
}

/* Borwein Thread para cálculo de Pi */
void *bot_calculate_pi(void *arg)
{
    bo_param *p = (bo_param *) arg;
    int j0, j1;

    /* "Bolhas" do método pipelined*/
    if(p->i < 2) return NULL;

    j0 = INDEX2(p->i-2, 0);
    j1 = INDEX2(p->i-2, 1);

    bo_pi_value(p->pi[j1], p->a[j1]);

    /* Se houve convergência, a variável stop é setada como true */
    if(has_converged(p->pi[j0], p->pi[j1])) p->stop = B_TRUE;

    return NULL;
}

